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Abstract 

The recently proposed Russian doll BCS model provides a simple example of a many body system 
whose renormalization group analysis reveals the existence of limit cycles in the running coupling 
constants of the model. The model was first studied using RG, mean field and numerical methods 
showing the Russian doll scaling of the spectrum, E n ~ Eq e _An , where A is the RG period. In 
this paper we use the recently discovered exact solution of this model to study the low energy 
spectrum. We find that, in addition to the standard quasiparticles, the electrons can bind into 
Cooper pairs that are different from those forming the condensate and with higher energy. These 
excited Cooper pairs can be described by a quantum number Q which appears in the Bethe ansatz 
equation and has a RG interpretation. 

PACS numbers: 74.20.Fg, 75.10.Jm, 71.10.Li, 73.21.La 
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I. INTRODUCTION 



The existence of limit cycles in the renormalization group (RG) is a possibility considered 
long ago by Wilson in the framework of High Energy Physics but only recently has it 
bund a concrete realization in several models in nuclear physics 0, quantum field theory 
HQ], quantum mechanics 0], superconductivity 0, S-matrix models 0,0], Bose- Einstein 
condensation [9( , effective low energy QCD [lOj , few body systems and Efimov states 
etc (for a review of see jllj].) The subject of duality cascades in supersymmetric gauge 
theory 2| is also suggestive of limit-cycle behavior. Chaotic flows have also been recently 
considered 0, [k| . The concepts of discrete scale invariance and quantum groups with q real 
are also closely related 

A RG limit cycle means that the coupling constants of the model are invariant under a 
finite RG transformation. There are generically two types of RG limit cycles: infrared and 
ultraviolet. In the former the limit cycles appear in the RG flow towards low energvand 
imply peculiar scaling properties of the spectrum termed as Russian doll scaling in p for 
obvious reasons. In particular if there are bound states, their energies E n (n — 0, 1, . . .) will 
scale as e~ nX where A is the "time" needed to complete an RG cycle. The ultraviolet RG 
limit cycles appear in the RG flow towards high energy and lead to log-periodic behavior 
of the scattering as a function of energy, and/or Russian doll behavior in the masses of 
resonances M n ~ e nX if they are present. 

Reference p proposed a slight modification of the BCS model of superconductivity, 
referred to here as the RD model, whose RG analysis revealed the existence of infrared limit 
cycles. The standard BCS model is given by the sum of a kinetic Hamiltonian describing 
the propagation of free electrons, plus a pairing Hamiltonian describing the scattering of a 
pair of electrons occupying time reversed states, say (k', T) an d ( — k', j), into another pair 
of states, say (k, f ) and (— k, |), with an amplitude Vk,k' lia [l6|. For s-wave pairing one 
can approximate the matrix element Vk,k' by a constant value G for all the incoming k' and 
outgoing momenta k within the Debye shell around the Fermi surface. To perform an RG 
analysis one defines a dimensionless BCS coupling constant g = GN(0) with N(0) the energy 
density at the Fermi level. Then under the RG flow g increases towards the IR becoming 
infinite at a scale signaling the formation of Cooper pairs. The modification added in p was 
to let the amplitude V^k' to pick up an extra imaginary piece proportional to ±ih where the 
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sign depends on wether the outgoing pair has greater or lower energy than the incoming pair 
(h is a dimensionless coupling constant). This of course breaks the time reversal symmetry 
which is a characteristic feature of the RD model. It was shown in [6j that the coupling h 
remains invariant under the RG flow while g exhibits a cyclic behavior with an RG period 
given by A = ir/h. Under this flow the coupling g becomes infinite at a finite scale but its 
value can be continued through minus infinity until it reaches its original value after the RG 
period A. One may expect from this behavior the existence of infinitely many scales related 
by A. Indeed, using the mean field BCS ansatz, which is valid for generic choices of the 
pair amplitudes Vk,k', one can show that the gap equation admits infinitely many solutions 
characterized by the gaps A n = A e~ nX (n = 0, 1, . . .), where the ground state corresponds 
to the solution with the highest value of the gap, A , and the other solutions, A n>0 , are high 
energy collective excited states. Let us call the latter solutions of the BCS gap equation 
the "superconduting dolls" or simply "dolls" and the integer n that label them the "nesting 
number" . The larger the nesting number, the larger is the size of the Cooper pairs forming 
the collective state, which is given by the correlation length £ n = £o e nA . 

These results showed in a simple case the intimate relationship between the cyclic proper- 
ties of the RG flows and the spectrum of the theory, which were also confirmed numerically 
in the case of one Cooper pair p. The case of more pairs was more difficult to deal with 
numerically. The DMRG or other ground state numerical methods were not helpful since 
one needs to explore high energy excited states to identify the "dolls". Fortunately the 
Russian doll BCS model has been recently shown to be exactly solvable by Dunning and 
Links using the Quantum Inverse Scattering Method [17]. This important result will enable 
us to explore in detail the spectrum of the RD model by solving the Bethe ansatz equations 
(BAE). This is the aim of this paper, whose organization is as follows. 

In section II we generalize the definition of the RD model, as given in [6], showing that 
the manifold of solutions of the gap equation with Russian doll scaling is a generic feature. 
In section III we review the Dunning and Links exact solution and show its agreement in 
the large N limit with the mean field solution obtained in section II (N is the number 
of particles). In section IV we focus on the solution of the BAE for the one Cooper pair 
problem, which serves to clarify some conceptual and technical issues before considering 
the many body case. In section V we explore numerically and analytically the low energy 
spectrum of the model in the large N limit, finding new types of elementary excitations 
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which are absent in the standard BCS model. These new kinds of elementary excitations 
can be characterized by a quantum number Q which has a clear meaning in both the RG 
and in the exact Bethe ansatz solution. In appendix A we review the analytic techniques 
used in section V and in appendix B we describe the elementary excitations of the standard 
BCS model in the canonical ensemble to facilitate the comparison with the excitations found 
in the RD model. 



II. THE RUSSIAN DOLL BCS MODEL: MEAN FIELD SOLUTION 



The model defined in |6| was based on a modification of the BCS model used to described 
the ultrasmall superconducting grains discovered by Ralph, Black and Tinkham H (see Q 
for a review). The latter model, also known as the picket-fence model in Nuclear Physics, 
has doubly degenerate discrete electronic energy levels which can be taken equally spaced 
or randomly distributed. It was claimed in p that the results obtained for the equally 
spaced model are valid in more general circumstances. We shall show in this section that 
this is indeed the case by considering a version of the RD model that is more similar to the 
standard formulation of the BCS model. 

Let us call c kya (c' k CT ) the destruction (creation) operator of an electron in the state 
with momenta k and spin a =t? I- The BCS pairing Hamiltonian (also called reduced 
Hamiltonian) is given by 

H =) e k n ki(T + 2J Vk,k' b[ by, (1) 

where e k is the energy of an electron with momenta k, = c\ a c k ^ is the number operator 
of the state (k, a), and b\ and b k are the creation and destruction operators of a pair of 
electrons in the states (k, ]) and (— k, {) 

b k = c I,T c -fc,|' hk = c ~Kl c k,V (2) 
The BCS ansatz for the ground state in the grand canonical ensemble is given by 

i^o>=n(^+^ & i)i°>> ( 3 ) 

k 

where |0) is the Fock vacuum of the fermion operators and u k , v k are variational parameters 
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whose values are given by, 



k 2 V E k r k 2 \ E k 



E k = \JZl + \A k \\ = e k - fi + V ktk . (4) 

The gap function A^ and the chemical potential /i are found by solving the gap equation 

A k = -J2v k , k >^, A k = \A k \e^, (5) 



1E k i 
k> K 



and the chemical potential equation, 



^ = 2EKI 2 = E( 1 -#)' ( 6 ) 



k k v ft 

where iV e is the number of electrons. If V kjk > is real, so is the gap function A*., up to an 
overall phase factor which can be chosen equal to 1 (i.e. <p k = 0). In the cases where the 
system possesses particle-hole symmetry around the Fermi surface the chemical potential 
coincides with the Fermi energy ep- Measuring all energies relative to ep one can take [i = 
We shall suppose that this is the case. 
The Russian doll BCS model is defined in terms of the scattering potential: 

J G + irj sign(e fc - e w ) \e k \ and \t k >\ < u c 

— Vk,k' — \ {<) 

I otherwise 

where sign(x) is the sign function. This potential describes for G > the attraction of 
electrons within a shell of width 2uj c centered around the Fermi surface plus an imaginary 
term which depends on the sign of the difference between the energies of the incoming and 
outgoing electrons. Since V k * k , = V k > )k the Hamiltonian is hermitian but the imaginary term 
breaks the time reversal symmetry. Setting rj = we recover the standard model used to 
describe s— wave superconductors Obviously the r\ interaction can be generalized to 

other type of symmetries such as p-wave, ci-wave, etc. Here we shall only consider the s-wave 
case. 

Let us next solve the gap equation (jSJ) for a large number of electrons where the discrete 
energy levels become a continuum. If we denote by N(e) the density of levels per energy 
then eq.(j2J), for the potential ((7|), becomes 

A(e) = £ de> N(e') [G + z V sign(e - e')] (8) 



where E(e) = ^e 2 + |A(e)| 2 if we assume particle- hole symmetry (// = ep — 0). Notice that 
Afc = A(e) depends on the momenta k through its energy e = e(k). In the standard BCS 
model where r\ — 0, eq.(|HJ) implies a constant gap, i.e. A(e) = A , whose value is given by 
the solution of the integral, 

— = de , v ; 9 

where we have used again particle-hole symmetry namely, i.e. N(e) = N(—e). Approximat- 
ing the density N(e) by its value at the Fermi level, N(0), one can perform the integral (JHJ) 
obtaining the well known result: 

2uj c e- 1 ' GN ^ (10) 



u sinh(l/GW(0)) 
The last expression is valid in the weak coupling case GN(0) < 1/4. 

When n ^ the gap A(e) depends on the energy. Differentiating eq.© with respect to 
e we find, 

de E(e) 1 j 

Hence the modulus of A(e) remains constant, A, and its phase varies with e, i.e. 

aw .a^, m = ,m (12) 

Without loss of generality we can choose 0(0) = so that 

* w -"jf'*§8 (13) 

The value of A is found by imposing eq.flSJ) at the Fermi energy, e = 0. Using eq. (fT2J) 
one can trade the integral over e into an integral over the phase </>, 

r<t>c A Ay Q 

1 = / _r (G - in sign(0)) e** = 1 - cos C + - sin C (14) 
J-4>c 2 V V 

where C = 4>{uj c ) is the value of the phase of the order parameter at the boundary of the 

Debye shell. The equation satisfied by 4> c is then 

tan0 c = C = Arctan ( ) + ttQ (15) 



G \G, 
where Arctan is the principal determination and Q is an integer. Plugging (fTK|) into (fTHJ) 
(for e = lj c ) one gets, 

1 



V 



Arctan (%) + ttQI = H de N ^ (16) 
VG/ Jo Je 2 + A 2 Q 



FIG. 1: Pictorial representation of two solutions of the RD equations with Q = (ground state) 
and Q = 1 (high energy state). The size of the Cooper pairs for Q = 0, £0, is smaller that the size 
of the Q = 1 pairs, £1, since ea.(|18|l implies £1 = e A £o- 



which is the equation that fixes A = Aq as a function of Q, G and r/. The positiveness of 
the RHS of this equation implies that Q must be a non negative integer. For each value of 
Q = 0, 1, . . . eq. (|TT)jl is identical to eq. © with an effective value of the coupling constant 
G Q given by 



1 1 

Gq V 



Arctan(^) + vrQ , Q = 0, 1, . . . (17) 



Hence in the weak coupling regime the gaps Aq satisfy the Russian doll scaling 

A Q ~ 2u c e- l / G ^ ~ A e- A< 3 (18) 

where 

A = p h = V N(0) (19) 
n 

can be identified with the "time" it takes the RG to complete a cycle jfjf • In the limit 77 — ► 
we get that Gq=o — > G and the solution Q = coincides with the unique BCS solution 
while the states with Q > have a vanishing gap and approach the asymptotically the 
Fermi state. 

The ground state of the RD model is the solution with the lowest total energy or equiva- 
lently with the lowest condensation energy E c . The latter quantity is defined as the difference 
between the energy of the ground state and the energy of the Fermi state, and it is given 
for the BCS state in the weak coupling regime by E c = — |iV(0)A 2 . Hence the ground state 
corresponds to the Q = solution, while the states with Q > are high energy excited 
states (see fig. |TJ). 

The quantum number Q not only determines the modulus of the superconducting order 
parameter Aq but also its phase. Eq. (jl5j) implies that the phase variation of 0(e) from the 



bottom to the top of the Debye shell is given by 2irQ, up to a constant, so that Q is a 
sort of winding number. This interpretation will be confirmed by the exact solution. These 
results generalize those obtained in to arbitrary energy densities N(e) and therefore are 
valid in two and three dimensions. In the ID case considered in |6( N(0) is given by 1/d, 
hence g = GN(0) = G/d and h = r]N(0) = rj/d, are dimensionless couplings. In more 
general cases where the matrix element Vk,k' depends on the momenta we also expect a 
RD behaviour as long as time reversal is broken, V k * k , = Vk',ki the reason being that the 
Physics is dominated by the vicinity of the Fermi surface in which case the simplified model 
(jZJ) is a good approximation. This points towards the universality of the RD model, whose 
experimental realization seems a priori feasible. 



III. RELATION BETWEEN THE EXACT AND MEAN FIELD SOLUTIONS 

We shall next summarize the exact solution of the Russian doll BCS model obtained 
by Dunning and Links We will use the notation adapted to the study of ultrasmall 
superconducting grains, although the results are more general as we have shown in the 
previous section. For grains, the single particle basis of momenta and spin is not appropriate 
due to the lack of translational invariance. However one can still use time reversed states 
\j, ±) created and destroyed by fermion operators c] ± and Cj t ± where j = 1, . . . , N labels N 
discrete energy levels Ej. The energy Ej represents the energy of a pair of electrons in a given 
level (it is twice the single particle energy of the previous section). In the equally spaced 
model the distance between two consecutive levels is fixed to 2d, i.e. Sj + i — Ej = 2d, so that 
—uj < Sj < uj where uj = Nd is twice the Debye energy u c . Henceforth all the energies will 
be twice their standard values. 

Let bj = Cj-Cj i+ , bj = c| + ct _ denote the usual pair operators. Using this notation the 
Hamiltonian (JTJ becomes 

N N 

H = Y, £ i b fa + E v '-r ( 2 °) 

/ 1 J:i' I 

where Vjy is the scattering potential. Similarly to eq.(JZJ), the RD model is defined by: 

— Vjji = G + i-r] sign(£j — Sj'). (21) 

In (j21j) we are assuming that all the pair energies Ej are different and labeled in increasing 
order, i.e. Sj < Ej> if j < j f . Hence the second term in (}2"T]) is equivalent to sign(j — j') so 
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that (j20j) can be written as 

N N 
j=l j>k=l 

where 

a = Arctan (£) , G = a/G 2 + n 2 (23) 

i — i 

The hamiltonian (J22)) is basically the one studied in [17j, where it was shown to be exactly 
solvable using the Quantum Inverse Scattering Method (QISM). These authors showed that 
(|2*2*|) appears in the transfer matrix of a certain vertex model as the second order term in a 
series expansion in the inverse of the spectral parameter. 

The standard BCS model, which corresponds to the limit rj — > of the RD model, is 
also exactly solvable and integrable. Its exact solution was obtained long ago by Richardson 
221 1 and the integrals of motion found much later by Cambiaggio, Rivas and Sarraceno 
24j). These results were rederived in the framework of the QISM 
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23| (for a review see 

in terms of an inhomogeneous XXX vertex model with twisted boundary conditions 
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281 ] . In this approach the hamiltonian and the conserved quantities appear in a 
semiclassical expansion of the transfer matrix in a parameter, 77' — ^ 0, that enters in the 
XXX vertex i?(«)-rnatrix defining the model. Expanding the i?(u)-matrix in this parameter, 
i.e. R(u) = 1 + 7]' r(u) + 0(r]' 2 ), yields the classical r(w)-matrix which satisfies the classical 
Yang-Baxter equation, while the i?-matrix satisfies the quantum Yang-Baxter equation. As 
we shall see below the latter parameter 7/ can be identified with the parameter 77 of the 
Russian doll BCS model. Hence, from the viewpoint of quantum integrability, the RD 
model is the "quantum" version of the standard BCS model which arises in the semiclassical 
limit: 

lim RD model = BCS model (24) 



Let us summarize the exact solution of the hamiltonian (J2UI22)) obtained in 
eigenstates in the sector with M electron pairs are found by solving the BAE, 



17j. The 



2i.a TT E a ~ £ j + i V TT E a -E b + 2ir] 



3 — = TT a & ^ ' (25) 

11 E a -Si -in 11 E a -E b -2iri y ' 

j=l a 3 1 6=l(^a) ' 

where the "rapidities" E a (a = 1, . . . , M) give the total energy of the state as, 

M 

E = E * ( 26 ) 
0=1 



In the limit where n — > we see from ()23j) that a ~ rj/G and then (|2*5j) becomes 

1 ^ 1 ^ 2 

• i - L -'a c 7 , iris - L -'a - L -'o 

These are the Richardson equations whose solution give the eigenstates of the BCS Hamil- 
tonian (i.e. 77 = 0). To solve the BAE (|2*3jl one can proceed as in the study of the antiferro- 
magnetic Heisenberg spin 1/2 chain by first taking the logarithm of the equations. Choosing 
the branch of the logarithm, 

^ log = Arctan ( ~) (28) 

we obtain from J2 



2i x — irj \x 



a + 7rQ a + J2 Alct & n ( E V _ )~ E Arctan ^ ) = ( 29 ) 
j=i ^ a ^ 6=1(^0) ^ a b ' 



where Q a (a = 1, . . . , M) are a set of integers. In the case of spin chains there is also a 
set of integers Qj associated to the rapidity variables whose choice determines the ground 
state and the excitations. For example the ground state of the antiferromagnetic spin chains 
corresponds to increasing Q's {Qj — j + const.), while holes in that distribution correspond 



29). 



to the elementary excitations of the model, i.e. the spinons 

Let us first show that the mean field solution of the Russian doll BCS obtained in section 
II coincides with the exact solution to leading order in N. We first need to recall that the 
mean field solution of the BCS model agrees to leading order in iV with Richardson's exact 
solution. This was proved by Gaudin j^O] and Richardson J^jJ who solved equation ()27j) in 
the large iV limit using an electrostatic analogy (see [3^] for the comparison between the 
analytical and numerical solutions, which is reviewed in the Appendix A). In the Gaudin and 
Richardson approaches the eq. ([27]) yield the BCS gap equation Q. The thermodynamic 
limit amounts to letting the number of levels iV go to infinity and the energy spacing d go 
to zero while u = Nd is kept constant. Eq.([27|) implies that E a ~ NG ~ Nd = u. On the 
other hand 77 ~ d ~ 1/N. Hence to leading order in iV eq. (j29|) becomes, 

1 N M 

(a + T<W + E _ E = „ (30) 

' j=l a 3 6=l(^o) 

The next to leading corrections in 1/N of this formula can be computed using the results of 
(3^, but shall not considered here. Taking Q a = Q for a = 1, ... , M, then eq. ([3T?]) becomes 
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the Richardson eq.(j2Z|) with a Q dependent value of the coupling constant 



a + irQ) = - Arctan + nQ (31) 



— = - (a + ttQ ) = - 
Gq V V 

where we have used eq. (J23|) . Hence the exact solution, to leading order in N, is given 
by the mean field solution with an effective value of the BCS coupling constant given by 
eq.(|31|). On the other hand this equation coincides with (|17jl. which shows that the mean 
field parameter Q, which labels the solutions of the gap equation, can be identified with the 
Bethe numbers Q a = Q, common to all the roots E a , which classifies the solutions of the 
BAE. This explains the origin of the quantum number Q in the exact solution, but it also 
suggests another possibility, namely the existence of states where Q a depends on the root 
E a . If we are interested in the low energy spectrum of the ground state it is clear that we 
must focus on the cases where Q a = for almost all values of a except for a finite and small 
number. 



IV. EXACT SOLUTION OF THE ONE COOPER PAIR PROBLEM 

As for the usual BCS model, it is useful to first consider the one Cooper pair problem, 
which is what we do in this section. Remarkably, the Bethe ansatz equations can be derived 
from the RG. 



A. RG derivation of the BAE 

The BAE f!29|) for a single Cooper pair (i.e. M — 1) with energy E reduces to: 



a + nQ + Arctan ( -— ^ — j = 
j=l V 3/ 



The solutions are the eigenvalues E of the Schrodinger equation 



(32) 



3-1 N 

( £ ._ G _ E )^ j = ( G + iv )J2ipt + (G-iv) J = 1,...,A (33) 

i=i e=j+i 

where ipj is the wave function of the one-pair state, i.e. X^iVV&j |0)- ft i s possible to 
derive eq.(|32j) using the RG method of Glazek and Wilson, which consists in the Gauss 
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elimination of degrees of freedom in a quantum mechanical problem P|. Using eq. ()33|) one 
first eliminates the high energy component ipN in terms of the low energy ones, 

IV 1=1 

Replacing this eq. back into eq. ()33|) yields an eq. for the remaining components, 

j-l N-l 

(sj-Gn^-E) if> j = (G N ^+ir]) J2^ + ( G N-i-iv) ^ j = l,...,N-l (35) 

1=1 f.=j+l 

where Gat_i is related to G and r\ by the eq. 

G 2 N + rf 

Gn~i = Gn H — , Gn = G (36) 

En — (jn — tL 

Notice that r\ remains invariant under the discrete RG transformation. Defining the quan- 
tities 

ft 71 

tana n = — , tan/3j = — (37) 

Cr n Hi — Ej 

eq.(|nnj) becomes, 

tanaAr_i = tan(aAr + (3n) => chjv-i = o>n + /3jv (mod it) (38) 

where a n is the principal part of the arctan (a n = Arctan(r]/G n )), etc. After p-RG steps 
one gets, 

p-i 

a N - p = a N + Y^ Pn-s + n Qp> 0<p< N -1 (39) 

where Q p is an integer. In the p = N — 1 step the Schodinger eq.(|3*5|) reduces to (ei — 
G\ — E)ipi = 0. The component ipi cannot be zero since otherwise the wave function would 
vanish. This implies that (e± — G\ — E) = 0. Using eqs. (}3l)j) and ()37|) the latter condition 
is formally equivalent to the eqs. G$ = oc and «o = 0. With the latter formal definitions 
eq.(JHHJ) can be extended to p = N, 

N 

= a N + ^(3j + nQ N (40) 

3=1 

which coincides with the BAE (|32j) upon the identifications = a and Qn = Q. 

The term ttQ in eq. (J4(J|) now has a RG interpretation. For low energy bound states, i.e. 
E = £ < Ei, eq. (jSll) becomes 

Arctan ( — — ] = Arctan ( 77- ) - 52 Arctan ( J + ttQ p (41) 

\G N -pJ \G N J ^ \£ N -j-e J 
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which can be viewed as a RG equation for the couplings. A cyclic RG flow is obtained if 
Gn-p = Gn, which happens under the condition, 

7tQ p = V" Arctan ( J (42) 

Q p can be identified with the number of RG cycles after p RG steps. The total number of 
RG cycles upon integration of all the degrees of freedom is given by 

N 



n ( c = Qn 



— \^ Arctan ( - 



v 



17 ~1 \ £ j ~ £ o 

where [. . .] stands for the integer part. For the equally spaced model we shall take Sj 
d(2j - N - 1) (j = 0, 1, . . . , N). Hence in the large N limit eq.flHJ) yields, 



(43) 



n 



h 



C ^ 



log AT, h = rj/d, A>1 (44) 



This equation was derived in p by taking the continuum limit of the discrete RG eq. (j3fij) 
which reads, 

f s =\^ 2 + h 2 ), g = G/d (45) 

where the RG scale s is given by s = \og(N /N) and N is the initial size of the system. 
The solution of this eq. 



g(s) = /it an 



T + tan \T 



9o = 9(0) (46) 



shows that the effective coupling g(s) is a periodic function of the scale s with a period given 
by Ai = 2n/h. In the many body case the period of the RG is a half of the one pair result, 
namely A = ir/h and the eq. (jjjj) has to be replaced by p 

h 

n c ~-logiV, iV >> 1 (47) 

7T 

Eqs. (I44jl and (|47j) admit a very simple derivation. First notice that the size of the system 
after a RG cycle is given by e~ x N. Hence in nc cycles, with 1 ~ e~ Xn °N, the system is 
reduced to a site. But in each cycle a bound state, in the one Cooper pair problem, or 
a "doll", in the many body problem (i.e. solution of the gap eq.), is produced and thus 
nc must give the number of bound states or "dolls". Eq. (J47j) also implies that in order to 
have at least one "doll", i.e. no > 1, the coupling h must be bigger than a critical value 
which depends on the size, h c = Ti/\ogN. If iV is for example the Avogadro number then 
h c = 0.057, which means that macroscopic samples may give rise to "dolls" even for small 
values of the /i-coupling. 
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B. Numerical solution 



The BAE (J32|) in the limit 77 — > becomes 



1 N 1 

~ = Y — 

G ^ E- 



(48) 



where we have choosen Q = in order to get a finite expression. Eq. (j48|) is the well known 
Cooper equation for a single pair and it has N solutions giving all the eigenstates of the one 
pair Hamiltonian (see figure EJ). 
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FIG. 2: Graphical solution of the ea. ()48j) for a equally spaced model for a coupling G = 1 and N = 
10 levels. |l6| |. The energies E are given by the abscisses of the points • which are the intersection 
of the horizontal dotted line —1/G with the function J2f=i ^/(^ ~ £ j) with Sj = 2j — N — 1 (d = 1). 
For G > the lowest value of E is always below e\ (Cooper pair state). 

Coming back to the RD model, eq. f|32|l . we show in fig. [21 all the solutions for a particular 
choice of parameters. Every solution, Ei, is characterized by a given value of Qi. Table 1 
collects these values for several examples. For small values of h all the Q's are zero as in the 
usual Cooper problem (not displayed in table 1). For h = 1 there is a single bound state 
with Q = and two high energy states with Q = —1. For h = 3 and 5 there are two bound 
states with Q = and 1. Finally for h = 10 there are three bound states with Q = 0, 1, 2 
(case depicted in fig El)- As shown in table 1, the number of bound states with positive value 
of Q is equal to the number of RG cycles (J43|) . This fact was explained above using limit 
cycle RG arguments. 
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FIG. 3: Graphical solution of eq. (|32j) for an equally spaced model with N = 10 energy levels 
(sj = 2j — N — 1, d = 1) and couplings g = G = 1, h = rj = 10. Every curve plots the LHS of 
eq . (|32j) for the values Q = —3, —2, —1, 0, 1, 2 from bottom to top. The solutions are the intersection 
of these curves with the real axis (points denoted by •). 



h 




Qi 


Q 2 


Q 3 


Q 4 


Q 5 


Qe 


Qi 


Q 8 


Q 9 


Qio 


1 





0* 























-1 


-l 


3 


1 


0* 


l* 














-1 


-1 


-2 


-l 


5 


1 


0* 


l* 


1 








-1 


-1 


-2 


-2 


-l 


10 


2 


0* 


l* 


2* 


l 





-1 


-2 


-3 


-2 


-l 



Table 1.- Values of Qi associated to the N = 10 eigenstates Ei for a equally spaced model 
with g — 1, h — 1, 3, 5, 10 and Ej — 2j — N — 1 (d = 1). The states are ordered increasingly 
with the energy. The symbol Q* means that it is a bound state, namely a state whose 
energy Ei is smaller that the lowest single pair energy E\ = —9. The second column gives 
the number or RG cycles as computed with eq. 
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V. THE MANY BODY PROBLEM 

The results of the previous sections strongly suggest the existence of a new type of exci- 
tations above the ground state of the RD model carrying non vanishing values of the Bethe 
numbers Q a . In the thermodynamic limit the BAE equation (J29|) becomes a Richardson like 
eq.(JHUj) which depends explicitly on the Q's, namely 



where 1/Gq = a/rj. In the rest of this section we shall study numerically and analytically the 
solutions of (J4T?|) for the equally spaced model with energy levels £j = (2j — N — 1) (d = 1). 

A. Numerical solutions 

The ground state of the RD model is given by the choice Q a = 0, Va. In the large N 
limit and in the strong coupling regime (i.e. go = Go/d > g c = 1.13459), the set of roots 
E a condense into an open arc T in the complex energy plane whose end points are given by 
a = £q_ — iA and b = Eo + iA, where Eo is the chemical potential and A is the value of the 
gap [32]. If the coupling g is smaller than g c , there is a fraction of roots that are real while 
the other ones are complex and form the arc. 

Next we present several numerical solutions of eq.(|49|) for a system with N = 100 and 
200 energy levels at half filling M = N/2. The BCS coupling constant is fixed to go = 2, so 
that all the roots for the ground state form complex conjugated pairs, except eventually for 
a single root which will be real. 

One pair: Figure |U shows the solutions E a corresponding to three choices: i) {Q a = 
0}M =1 , ii) Q l = 1 and {Q a = 0}*£ 2 and iii) Q 1 = 2 and {Q a = 0}*£ 2 . For the choice i) the 
M roots E a form the largest arc located to the left, while for ii) and iii) the M — 1 roots 
with Q a = form a smaller arc, while the real root with Q\ = 1 or 2 lie on the real axis 
far apart from the arc. We have found more than one solution with Qi = 1 or 2 which 
must correspond to higher excited states of the one Copper pair problem, as shown in the 
previous section. In table 2 we collect the numerical and theoretical values of the real root 
Ei = £i for several values of the ratio Qi/h, as well as the excitation energy of the state. 




a = 1, . . . ,M 



(49) 
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FIG. 4: Plot of the roots {E a }^ =l of eq.gU) for the three cases described in the text and the values 
N = 100, M = 51, go = 2,h = 2. The roots with Q a = 0,Va form the "ground state" arc which 
shrinks when Q\ = 1 or 2 (the roots forming the two arcs are related by the links • — •). The real 
root £i for Q = 2 is close to the bottom of the energy band uj = —Nd = —100. All the energies in 
figs. 4-7 are given in units of d. 

The theoretical values will be obtained in the next subsection. 



Qi/h 


fnum 
SI 


tth 
SI 


tth' 
SI 


E 


pnum 

exc 


771th' 

cxc 


3/2 


-97.365 


-100.021 


-100.021 


-2694.635 


217.998 


216.405 


1 


-100.277 


-100.477 


-100.481 


-2693.203 


219.43 


216.618 


1/2 


-112.317 


-111.941 


-112.059 


-2686.975 


225.658 


222.225 


1/3 


-132.24 


-140.931 


-141.476 


-2697.788 


214.845 


238.416 



Table 2.- Numerical and theoretical values of the real root £1, total energy E and excita- 
tion energy E exc = E — E (E is the ground state energy) for several values of Qi/h. The 
rest of the parameters are fixed to iV = 100, M = 51, g = 2. The cases depicted in fig. 
namely h = 2, Qi = 1 and 2, correspond to the values Qi/h = 0.5 and 1 in bold face. 

Two pairs: Fig. O-left shows the solutions E a corresponding to two choices: i) {Q a = 
0}aLi an d ii) Qi = Q2 = 1 and {Q a = 0}^f =3 In the second case the roots E^ 2 = ^1,2 form a 
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FIG. 5: Left: Roots E a of eq.(jl| for the values N = 100, M = 50, g = 2,h = 2. There are two 
roots £12 with = 1 forming a complex conjugated pair. Right: Different choices of the roots 
for which Q a = 1 yield the same result. 

complex conjugated pair near the real axis. Observe how the M roots of the ground state arc 
reorganize themselves into a new arc with two less roots having Qi^ = 1. In the numerical 
program one can choose another pair of roots, say E 7 and E$, getting the same result. This 
is shown in Fig El-right. In table 3 we collect the numerical and theoretical values for the 
complex roots E\^ = £ lj2 for several values of the ratio Qi^/h, as well as the excitation 
energy of the state. The theoretical values will be obtained in the next subsection. 



Qi,2/h 


tnum 
Sl,2 


SI, 2 


E 


rpmim 
^exc 




0.83 


-100.29 ± 1.515 i 


-100.092 ± 1.339 i 


-2481.341 


432.458 


432.869 


0.50 


-111.008 ± 5.058 i 


-110.632 ± 4.963 i 


-2470.368 


443.431 


442.935 


0.38 


-126.262 ± 11.498 % 


-126.457 ± 8.12 i 


-2464.325 


449.474 


459.445 



Table 3.- Numerical and theoretical values of the complex roots E\^ = £1^, total energy E 
and excitation energy E exc for several ratios Qi^/h. The rest of the parameters are fixed to 
iV = 100, M = 50, g = 2. The case depicted in fig. namely h = 2, Q 12 = 1, corresponds 
to the value Qi^/h = 0.5 in bold face. 
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FIG. 6: Roots E a of eq.jHJ) for the cases: Left: N = 200, M = 100, g = 2,h = 2. There are two 
complex roots with Q\ t 2 = 1 (the closest to the arcs) and a real root with Q% = 2 (the farther from 
the arcs). Right: N = 200, M = 101, g = 2, h = 4. There are five roots with Q a = 1. 

Three pairs: Fig.^J-left shows a solution with three roots formed by a complex pair with 
Qi2 = 1 and a real pair with Q% = 2. In table 4 we collect the results. 



fnum 
Sl,2 


fnum 
S3 


d-th 
SI, 2 


d-th 
SI, 2 


pnum 

exc 


E th 

exc 


-218.229 ± 7.365i 


-201.356 


-222.578 ± 7.115z 


-200.949 


1319.06 


1320.5 



Table 4.- Numerical and theoretical values associated to fig. El-left. 

Five pairs: Fig. El-right shows a solution with five roots with <5i,...,5 = 1, one of which is 
real and the remaining four are two complex conjugated pairs. As can be seen they form a 
small arc which one expects to become larger with the number of Q a = 1 roots. 

M/2 pairs: Fig. [7| shows a solution with M = 100 roots where half of them have Q a = 
and the other half have Q a = 1. This represents a high energy state which is intermediate 
between the Q = and Q = 1 states where all the roots condense into single arcs. 

The main conclusions we can draw from the previous numerical results, and others not 
shown above, are the following. In the thermodynamic limit where the number of roots E a 
with Q a > is kept finite as iV and M become very large we have: 
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FIG. 7: Roots E a of eq.® for the values N = 200, M = 100, g = 2, h = 4. The wider arc on 
the LHS is the usual ground state arc. The Y shaped arc in the middle corresponds to M/2 = 50 
roots with Q a = 1, while the two disconnected branches above and below it contain the remaining 
50 roots with Q a = 0. 

• The roots with Q a = form an open arc T which is a slight perturbation of the arc T 
formed by the M roots of the ground state. 

• The roots with a non vanishing value of Q a fall into arcs characterized by Q a . 



B. Analytic solution 

The previous results suggest the steps to follow in the analytical study of eq. (j4T?j) . We 
shall employ the methods of complex analysis developed to establish the connection between 
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the Richardson eqs. in the thermodynamic limit and the mean field BCS solution 
(se appendix A for a review). 

Let us first rewrite eq.tjUU) as follows, 

- 1/2 - 1 1 

E-hr- E .e=e-w = ° (50) 



3=1 3 b=l(^a) " 



where 



7T = 7T + > 7T = -Arctan - 51 

G a G V G r] \GJ 
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and let us suppose that Q a = for all roots E a except for a finite number in the limit 
N — > oo. We shall call the latter roots £ a (a = 1, . . . , D). Making an electrostatic analogy, 
eqs.([5L)j) and (|51j) imply that on each root E a acts an electric field —l/2G a , which depends 
on the value of Q a . The roots with Q a = see an electric field —1/2Gq, while the Q a > 
roots see a stronger field. Assuming that all the roots Q a = form a single arc T, then the 
roots £ a , with Q a > 0, must lie outside T, namely 

E.J^l Q - = ° (52) 

[tat r, g a > o 

There may also exist roots with Q a = lying outside the arc T. They will be considered 
later on. In the continuum limit eqs.([50|) split into two sets of eqs. (see appendix A for 
definitions), 



d£p(c) t Ar-pf.^-^ e.*f (54, 



where the first equation holds for the roots £ on the arc T and the second equation holds 
for the roots £ Q outside T. The density of roots r(£) along the arc T, and the total energy 
of the state E are given by (see eqs. (|8H82j) ) 



j \d£\ r(0 = M-D, (55) 

D 

/ K|£r(0 + ^^ = ^- 

Jr a=l 



(56) 

The solution of eqs. (|53l54)) follows the same steps as in Appendix A. Let us summarize 
the results. The density function r(£) can be found from a formula similar to (JSlj) with a 
modified electric field h(£) given by, 

with 
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R(Z) = y/it-aKt-b) (58) 

and where a = So — iA, b = Eq + iA are the end points of the arc V . Using eq. (|57jl into (|53|) 
one gets after some algebra, 

1 _ f depie) 1 _ m 



2G 



o Jn 



which is a modified gap equation analogous to ()91j). Similarly (|55|) gives the modified 
chemical potential eq. (see (jSHJ)), 



M-D = /&p(s) (1- , )-E(l- , 1 (60) 

( e _e )a + AV °=i V - £ ) 2 + A* 



These two equations determine the end points £o ± of the arc V . They are formally 
equivalent to the eqs. (}9Tj) and (|9*3j) if one defines an effective density 



D 



p(e) = p{e) + 5p(s), 5p{e) = - £ *( e - £ a ) (61) 



a=l 



Similarly the eq. (J53|) for the roots £ a turns into, 



D 



^ = [ d£ A^LSM - T 1 ggg) , Rfoo m 

2v Jn e-tiaR(e) & ~ ^ Rfo) R{£ a ) 

where the RHS can be identified with the electric field h(£ a ) (see eq. ()57|l ). after subtracting 
the pole at £ = £ a . Finally, the energy of the state can be derived from eq. lJHfi)) and it reads, 



E = ~+ [dee 1- - £ £o =\p(e) (63) 



4G 



n 



D D 



e - e y + 



a=l \ ^(^-e )2 + A2 

Doing a computation similar to the one that leads to eqs. (jlOlj) and ()103|) in Appendix B 
one gets the excitation energy of the state E exc = E — E (E is the ground state energy) 
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D 



E exc = J2 \/(6>-£o) 2 + A 2 , (64) 



0=1 



where A and Eo are the unperturbed values of the gap and chemical potentials given by the 
solution of eqs.(jnU) and (|9"5|) with G = G . 

C. Analytic versus numerics 

Let us compare the analytic and the numerical results obtained previously. For one real 
pair £1 the eq.(jB2J) becomes, 

^=fde^m + m (65) 
2?7 Jn £ - 6 R(e) R(^) 

In the large N limit we can approximate R(£) by R{£), so that (|u3jl simplifies, 

W = K(l) + ({.-*.)' + *» (66) 
where we have used ([93)1 . In the equally spaced model at half filling the integral giving h(£) 
is given by the formula j^ . 

1 ({-,)(^-^ +v( ^ + a(A- + ^)) 

4d (£ + w)(A 2 + ?w + V(A 2 + £ 2 )(A 2 + w 2 )) 

To compare with the numerical results collected in table 2 we solve eq. ()66|) for the values, 
N = 100, g = 2, d = 1 which implies that A = Ad/ sinh( 1/5-0) = 191.903, e = and 
u; = dN = 100. The result is given in table 2 in the column £j; h '. The excitation energy is 



obtained using (EH) , namely £ exc = v / ^FF + A 2 and the values 

appear in the column £^ 
of table 2. Observe the good agreement between the numerical and theoretical values. 

As can be seen from eq. (|67J) the field h(£i) is proportional to 1/d ~ N, hence in the large 
N limit the eq. (JHHJ) can be approximated by 

^ = Kb) (68) 
On the other hand the change of variables j^] 

i A cosh u , , 

* = / . 2 (69) 

1 — ( — sinh-u) 
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leads to a very simple expression for h(£), namely 



So that the solution of eq. (}68|) is given simply by 

Acosh(7rQ//i) 
U + (£sinh(7rQ//>)) 2 

The column £j h of table 2 collects several values of (fTTj) . which are quite close to the values 
of £* h '. The solution (JUJ) remains real and below the bottom of the band, i.e. £1 < — lu, 
whenever irQ/h > 1/go. 

To check the results of table 3 one has to solve eq. ()62|) for two complex roots, say £1 and 
£ 2 = 61 with Q 1 = Q 2 , 

^ = %) > §4, fa^ ({2) -^f4 (72) 

277 £ 2 - 6 #(60 277 & - £ 2 R(€i) 

Notice that we drop the term R'(£)/ R(£) in these equations. The numerical results are shown 
in the column £j h 2 . Similarly the excitation energy is given by E exc = 2Re[\J (£j h ) 2 + A 2 ] 
Finally, the theoretical results of table 4 are obtained by choosing £j h 2 as the solutions of 
eqs.([72^1. and ^ h from eq.(J7TJ). Strictely speaking one must solve eqs.ljo^j) for three roots 
with Q a = 1,1,2, but they can be approximated as above. The energy is given by 
summing over the three roots £* h . 



D. The Q-excitations in the PBCS ansatz 

Let us summarize the results obtained so far. Using the grand canonical BCS ansatz we 
showed the existence of different solutions of the gap equation in the RD model characterized 
by an integer Q > 0, where the Q = solution corresponds to the ground state and the 
solutions with Q > 1 correspond to higher excited states with different condensation energies. 
Later on, the exact solution of the model led us to identify the integer Q with the Bethe 
numbers appearing in the solution of the BAE for the M roots, namely Q = Q a , Va. This 
suggested the existence of low lying excited states for which Q a = for most of the roots, in 
the thermodynamic limit, except for a finite number of them, for which Q a > 1. These new 
type of excitations cannot be derived from the standard mean field analysis in the grand 
canonical ensemble which are given in terms of the familiar Bogoliubov quasiparticles. This 
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is why we have to use the exact solution to find them. One may wonder however if there 
is an alternative way to derive the Q— excitations, which could be valid in more general 
circumstances. We now argue that this is indeed possible using the projected BCS ansatz. 
Let us write the BCS state defined in eq.(J3J) as 

|Vo)cxexp^(fo&n |0> (73) 

where = Vk/u^ is the wave function of the Cooper pair in momentum space. A state with 
M pairs can be obtained from ()73|) keeping the M th -term of the Taylor expansion, i.e. 

|M)<xfe>6^ |0> (74) 

which is called the Projected BCS state (PBCS). While the grand canonical BCS state (J73J) 
is an eigenstate of the phase operator we see that (J74"Jl is an eigenstate of the electron number 
operator. In the large N limit the results obtained with both ansatze agree to leading order 
in N. However for finite values of N the BCS and PBCS give different results as ocurred 
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when they were applied to the study of ultrasmall superconducting grains 
the RD model the PBCS state corresponding to the <5 th -solution can be written as 



35, 361. In 



/ \ M Q 

|M Q >ocr>>f&n |0> (75) 

where is the associated wave function. The ground state is given by (J75J) with Q — 0. We 
shall conjecture that the states found in the subsections V-A and V-B can be approximated 
by the following PBCS states 

t\ M o/ \ M i 

X>i 0)& i) (E^ 6 *) ••• i°) (76) 

where Mo ~ O(M) is the number of roots with Q a = 0, M\ ~ 0(1) the number of roots 
with Q a = 1, etc. Eq.()76J) certainly holds in the case where the operators bk are ordinary 
bosons instead of hard core ones. In this case the RD Hamiltonian can be diagonalized by 
a canonical transformation of the boson operators which amounts to solving the BAE for 
one pair, eq.([32|). Each term in (J76j) would correspond to the bound state solutions of the 
one pair problem with Q = 0, 1, . . .. 

An important consequence of this discussion is that the Q-excitations must behave as 
ordinary bosons. In the canonical picture one simply adds one more pair to the Q > 
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states, which then increases the excitation energy. In the exact solution this corresponds to 
removing pairs from the ground state arc and place them into the small arcs with Q a > 0. 
This property is in sharp contrast with the standard quasiparticles which are fermions. 



E. Pair-hole excitations as Q a = states 

So far we have considered the solutions of eq. (J49j) where some of the Q a are non zero 
which has the effect that the associated roots come out from the ground state arc. In 
it was shown that the Richardson eqs. (J27|) already contains many solutions where this 
happens. The latter solutions correspond to the pair-hole excitations of the BCS model, to 
be distinguished from the excitations where the energy levels are blocked by single electrons, 
(see appendix B) . In the large N limit the pair- hole excitations are characterized by the fact 
that Ng roots E a lie outside the arc Y formed by the M — Ng remaining ones j^. m figure 
IH1 we show an example of two excited states with Ng = 1 and 2 roots for a model with 
M = 20 roots. Figure is similar to figlU since in both cases one root comes out of the 
ground state arc T. They differ however in the value of Qx, which is 1 in fig. H]and in fig. 



IHK- In [37[ it was found that the roots E a (a — 1, . . . , Ng), lying outside the arc T satisfy 
in the large N limit the following eq. 

n= /V p(g) ggo) y i R&) | (77] 
Jn e-i a r {£ ) f£ (fr - Q Rfo) R(£*) ' 1 ) 

which coincides with eq. 1)62)1 upon the choice Q a = 0. Moreover the energy of these excited 
states is given by, 

N a 

E exc = V^-^ + A 2 , (78) 

a=l 

This formula coincides with eq. ()64jl giving the energy of the Q Q -excitations. Hence the 

rn 

pair-hole excitations of [37] can be seen as Q a = excitations. 

Summarizing, in the RD model there are three types of excitations in the canonical 
ensemble, i) the ones obtained by the blocking of energy levels, ii) the pair-hole excitations 
and iii) the Q— excitations. The first two have already been considered for the standard BCS 
model and in the g.c. ensemble they correspond to the quasiparticles which have fermionic 
statistics. The last ones are specific of the RD model and have bosonic statistics. These 
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FIG. 8: Position M = 20 pairs of the states at g = 1.5. The arcs T for 19 pairs (left) and 18 pairs 
(right) are a slight modification of the ground state arc T (20 pairs) 



results must have profound consequences in the thermodynamical properties of the system 
and its response to external fields. 



VI. CONCLUSIONS AND PROSPECTS 



In this paper we have analyzed in more depth the Russian doll BCS model using its exact 
solution. We have shown that the mean field solution obtained in [f| agrees in the limit of 
large number of particles N with the exact solution to leading order in N. In doing so we 
have identified the integer Q, that labels the mean field solutions, with the Bethe numbers 
appearing in the BAE. This integer also has a RG meaning since it counts the number of 
RG cycles. This idea is clarified by deriving the BAE for one Cooper pair using the RG 
method of Glazek and Wilson 5]. The numerical solution of the latter equation shows the 
appearance of bound states with Q > and of unbounded states with Q < 0. We have also 
studied the many body case in the large N limit, where the BAE can be approximated by 
Richardson like equations where the value of the effective BCS coupling constant depends 
on the Bethe numbers Q a . We have solved numerically and analytically these equations 
showing the existence of solutions where the Bethe numbers Q a may vary with the roots. 
The most interesting case is when Q a vanishes for most of the roots except for a small number 
where they are positive. These solutions correspond to low energy excitations where Cooper 
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pairs, forming the ground state, jump into excited states characterized by positive values of 
Q a > 0. In this sense Q a becomes a sort of principal quantum number of the Cooper pairs. 
These new type of elementary excitations should be added to the standard ones to describe 
the low energy spectrum of the model, whose thermodynamical properties and response to 
external fields will be modified. 
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Appendix A: Gaudin's electrostatic solution of Richardson equations 

In this section we review the electrostatic method employed by Gaudin in the proof of 
the asymptotic agreement of the exact and mean field solutions (see also 2, 38j|). The 



starting point of Gaudin's approach is the observation that the Richardson eq.(|2*7|). that we 
shall rewrite as 

^ 1/2 ^ _J 1_ 

can be seen as the equilibrium conditions for a set of mobile charges +1 located at the 
positions E a in the complex plane under the effect of a constant electric field — 1/2G and 
another set of charges —1/2 at the positions Ej. In the large N limit the pair energy levels 
Ej will be equivalent to a negative charge density —p{s) located on an interval Q = (—u,u) 
of the real axis. The total charge of this interval is given by 

N 

r at 1 

- / de p{e) = -- , p(e) = - £ 5(e - e s ) (80) 
Jn j=1 

The numerical solutions of eqs. ()79|) are either real or complex roots in which case they 
come in conjugated pairs 0, 13 m the large N limit the complex roots form an open arc 
T whose end points are given by a = Eq — iA and b = £q + iA, where eq is the chemical 
potential and A is the gap (see fig. EJ). Let us assume for simplicity that all the roots are 
complex. We shall call r(£) the linear charge density of roots E a along the arc V. Hence, 
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FIG. 9: Plot of the roots En for the equally spaced model in the complex plane (taken from 
The discrete symbols denote the numerical values for M = 100. The continuous lines are the 
analytical curves obtained in 32]. All the energies are in units of oj. 



the total number of pairs, M, and energy, E, are given by 



m r(0 = M, 



The continuum limit of eqs. ()79|) is 



de p(e) 



P 



2^ = °' 



(81) 
(82) 

(83) 



which implies the vanishing of the total electric field on every point of the arc T. The 
solution of eqs. (|H3*j) can be found as follows. First of all, let us orient the arc T from the 
point a to the point b, and call L an anticlockwise path encircling T. We look for an analytic 
field h(£) outside F and the set Q, such that 

1 



r(0 m = ^(MO - MO) ^r, 



(84) 



where h + (£) and h-(£) denote the limit values of h(£) to the right and left of T. This can 
be understood using the electrostatic equivalence, considering that the electric field presents 
a discontinuity proportional to the linear charge density when crossing the arc T. Next, we 
define a function R(£), with cuts along the curve T, by the equation 



(85) 
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and look for a solution which vanishes at the boundary points of T, in the form 

ho = m I ^ . (86) 

This field has to be constant at infinity, as can be verified explicitly. The contour integral 
surrounding the charge density r(£) in (jHHj) can be expressed as 

ki r(o i i df , h + (c) - h-(o = ±r^m (e cr (87) 



where Cr is the region outside the curve T. Using eqs. (J85j) and one finds for the 

principal value of (JHTj) 

(88) 

where we have deformed the contour of integration L into two contours, one encircling the 
interval Q (first term on the RHS) and another one around the infinity (second term). We 
are assuming that T, and consequently L, do not cut the interval Q, which happens int the 



equally spaced model when g > g c = 1.13459 
a solution is obtained provided 



32j |. Plugging eq. (jHHj) into (JEHJ), we see that 



Using the eq.(j2DJ) becomes, 



"wrh- (90) 



dep(e) _ 1 (91) 



'n - e„r- + A 2 2G 

which is nothing but the BCS gap eq.® in the appropriate normalizations (see eq.([97|) 
below). The field h(£) gives the charge density r(£), 



r(0 = -\h(0\, eer, (92) 

7T 



and its value is given by replacing (j89j) into (Joojl . i.e 



n 
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The equation fixing the arc T are the equipotential curves of the total distribution, 

K I d? h(g) =0, (er. (94) 

J a 

Similarly, eq. (|5T|) becomes the chemical potential equation 

M=^-ldt h{i) = [ de p(e) [ 1 - . £ = £ ° = ) , 
while eq. (jH2J) gives the BCS expression for the ground state energy, 



(95) 



E = i I ^ * MO = ~+ fdee [1- - n " f° ao 1 (90) 



Comparing these equations with the corresponding ones in the BCS theory, we deduce the 
following relations between A, e , and p(e), and A BC s (BCS gap), fi (chemical potential), 
and n(e) (single particle energy density): 

A = 2A BCS , e = 2p, p (e) = in(|). (97) 

Appendix B: Elementary excitations of BCS in the canonical ensemble 

In the grand canonical ensemble the excited states can be obtained acting on the GS 
ansatz |^ ) with the Bogoliubov operators 7j jCr (er = ±) 

lh,cn ■ ■ ■ 7i„,a„ I'M, 7i,± = u j c i± T Vj c] T , (98) 

where the variational parameters Uj, Vj are given by eqs. similar to (j3J). The excitation 
energy is given by (recall the factor of 2 in our conventions as compare to the standard 
ones) . 

E exc = E-E = ^J2 \j ( £ i ~ £ o) 2 + A 2 (99) 

3 

We have to distinguish between two sorts of excitations: El) quasiparticles occupying differ- 
ent energy levels, which means that all the j's in eq. (|9*Hj) are different, and E2) quasiparticles 
occupying the same energy level, i.e. 7j + 7j_. In the latter case the factor (uj + Vj &t) in the 
ground state is replaced by the factor (vj — Uj The states E2 are called "real pairs" to 
be distinguished from the "virtual pairs" that build up the ground state 4§] 

In the canonical ensemble the Bogoliubov operators do not even make sense since they 
change the particle number which is fixed by definition. It is however possible to see the 
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correspondence of these excitations in the canonical ensemble. The excitations of type El 
correspond to blocking of energy levels and we discuss them next in detail. The excitation 
of type E2 correspond to pair-hole excitations and they are discussed in subsection V-E. 
Blocking of energy levels 

Let us break P pairs of the M pairs forming the ground state and place the corresponding 
2P electrons in different 2P levels belonging to the set B = {ji < 32 < ■ ■ ■ < 32p}- All these 
levels are singly occupied and therefore the corresponding electrons decouple from the rest 
of the system contributing only with their free kinetic energy \£j, j G B. The remaining 
M — P pairs are then allowed to occupy the N — 2P energy levels that are left (see fig. El 
for an example). The problem we are left with is entirely similar to the original one after 




FIG. 10: Excited state obtained by breaking a "Cooper pair" at G = 0. On the RHS the singly 
occupied levels are blocked and decouple from the rest of the system. 

removing the blocked levels. This is turn is equivalent to consider the modified density of 
levels, 

p(s) = p{e) + 5p(e), 5p(e) = -5 J>( e " e i) ( 10 °) 

jeB 

where p(e) was defined in (|80)1. This small perturbation of p(e) translates into a small 
perturbation of the density, i.e. r(£) = r(£) + 5r(£) and that of the arc T whose ends points 
will move slightly, i.e. £q — £q + 5eq and A = A + SA. The eqs. for the deviations in the 
chemical potential and the gap can be obtained by taking the variation of the gap eq. (J91|) 
and the chemical potential eq.(P?5|): 

Jje^[A5A- { e-e )5e o] = jy-M (101) 

/ d£ 4rTT t( £ - £ o)5A + Afe ] = 5M- [ ds5p (l " : " 
Jn R{£r Jn V 



R{e) 
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where R(e) is given by eq.()85|). Taking the variation in ()96|) one finds 

5E = + / dee [(e - e )5A + A6e ] (102) 

+ / n d " W£) ( 1_ if) 

Using the gap eq. for the first term in the RHS and applying the eqs. (|1U1|) one gets, 

5E = e 5M + I de 5p{e) [e - e - R(e)] (103) 
Jn 

= 5 

where we used (|1UU|) and <5M = —P. We have to add to eq. (|103|) the contribution of the 
single occupied levels, namely \ J2jeB £ j which then coincides with the Bogoliubov formula 
(jnnj)- This proves the one to one correspondence between the Bogoliubov quasiparticles 
occupying different energy levels and the blocked levels in the canonical ensemble. 
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